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SPHRAY: A Smoothed Particle Hydrodynamics Ray 
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ABSTRACT 

We introduce the publically available code SPHRAY, a Smoothed Particle Hydrody- 
namics (SPH) ray tracer designed to solve the 3D, time dependent, radiative transfer 
(RT) equation for cosmological density fields. The SPH nature of SPHRAY makes the in- 
corporation of separate hydrodynamics and gravity solvers very natural. SPHRAY relies 
on a Monte Carlo (MC) ray tracing scheme that does not interpolate the SPH particles 
onto a grid but instead integrates directly through the SPH kernels. Given an arbitrary 
(series of) SPH density field(s) and a description of the sources of ionizing radiation, 
the code will calculate the non-equilibrium ionization and temperature state of Hydro- 
gen (HI, HII) and Helium (Hel, Hell, Helll). The sources of radiation can include point 
like objects, diffuse recombination radiation, and a background field from outside the 
computational volume. The MC ray tracing implementation allows for the quick intro- 
duction of new physics and is parallclization friendly. A quick Axis Aligned Bounding 
Box (AABB) test taken from computer graphics applications allows for the accel- 
eration of the raytracing component. We present the algorithms used in SPHRAY and 
verify the code by performing the test problems detailed in the recent Radiative Trans- 
fer Comparison Project of Iliev et. al. The source code for SPHRAY and example SPH 
density fields are made available on a companion website (www.sphray.org). 

Key words: cosmology, theory, numerical methods, N-body, SPH, ray tracing, Monte 
Carlo, simulations, radiative transfer, reionization, Stromgren 
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1 INTRODUCTION 



In numerical cosmology, prescriptions for the treatment of 
gravity and hydrodynamics are well developed and have 
been validated against one another in several comparison 
studies (see Frenk et al, 1999; O'Shea et al., 2005; Heit- 
mann et al, 2005, 2007; Regan et al, 2007; Agertz et al., 
2007; Price, 2007). The density and temperature fields they 
produce provide input for sub-resolution models of star for- 
mation and feedback via supernovae (e.g. Springel & Hern- 
quist, 2003) and black holes (e.g. Di Matteo et al., 2007). 
Numerical radiative transfer (RT) techniques, necessary to 
calculate the interaction of the ionizing photons produced 
by these sources with the cosmological gas, have not yet 
reached the level of maturity attained by A^-body and gas 
dynamics solvers. Flexible and accurate RT techniques, vali- 
dated against analytic solutions and in comparison projects, 
are necessary to properly interpret many observations and 
guide the development of theoretical models from cosmolog- 
ical through stellar scales. This is especially true for analy- 



sis of upcoming 21 cm surveys such as 21CMA 1 (formerly 
PAST), LOFAR 2 , MWA 3 , SKA 4 ; modeling absorption 
lines in the spectra of high redshift quasars and gamma ray 
burst afterglows, and understanding the feedback processes 
which influence star and galaxy formation. 

The introduction of 3D radiative transfer into cosmo- 
logical simulations is complicated by several issues. The spe- 
cific intensity I v = /(x, n, v, t) is a function of seven vari- 
ables leading to a solution space with high dimensionality. 
R-type ionization fronts can travel at nearly the speed of 
light through underdense regions and many times the speed 
of sound in dense regions leading to radiative time scales 
orders of magnitude smaller than dynamical time scales. In 
addition, radiative transfer and hydrodynamic processes are 
coupled. For example, photo heating creates large pressure 
gradients near luminous sources and can modify star forma- 
tion rates while hydrodynamic temperature changes affect 

1 http: / /21cma. bao.ac.cn/index.php 

2 www.lofar.org 

3 www.haystack.mit.edu/ast/arrays/mwa 

4 www.skatelescope.org 
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the recombination, collisional ionization, and radiative cool- 
ing rates of the cosmological gas. 

In designing a numerical radiative transfer scheme, it 
is practical to utilize the hydrodynamic frameworks that 
have already been developed. These can generally be catego- 
rized as Lagrangian, particle based methods (see Monaghan, 
1992, for a review of Smoothed Particle Hydrodynamics) or 
Eulerian, grid based methods (see e.g. Norman, 2004, for 
information on Adaptive Mesh Refinement codes). In this 
paper, we describe SPHRAY , a code that performs radiative 
transfer calculations on Smoothed Particle Hydrodynamics 
(SPH) density fields. A flexible, general purpose, radiative 
transfer method tightly coupled to SPH hydrodynamic sim- 
ulations could be adapted to handle many different astro- 
physical problems. Currently, SPHRAY works on static den- 
sity fields. However, it calculates quantities that are equiv- 
alent to the change in specific energy (or entropy) for indi- 
vidual SPH particles due to photoionization/photoheating. 
This makes the coupling of SPHRAY with gravity and hydro- 
dynamics relatively straightforward. We leave this to future 
work. 

SPH (Lucy, 1977; Gingold & Monaghan, 1977) is a grid- 
less, Lagrangian method that discretizes a fluid into parti- 
cles. The self-gravity of these particles can be treated in the 
same way as Af-Body particles, but they are also subject 
to hydrodynamic forces. The combination of SPH with tree 
structures (Hernquist & Katz, 1989) and of tree structures 
with the particle-mesh method (Xu, 1995; Bagla, 2002), pro- 
vides a very flexible computational tool. 

The earliest combinations of SPH and RT were by Lucy 
(1977) - one of the papers that introduced SPH - and Brook- 
shaw (1985). These authors modeled radiation transport as 
a diffusion process. The study of dense astrophysical regions 
such as molecular clouds and collapsing protostars has con- 
tinued along this line in the work of Whitehouse & Bate 
(2004); Whitehouse et al. (2005); Viau et al. (2006) and 
Mayer et al. (2007). 

The highly variable optical depths through voids and 
Lyman Limit systems in cosmological volumes do not lend 
themselves to a diffusion description of radiation. Here, ei- 
ther raytracing schemes or moment methods must be used. 
SPHRAY uses Monte Carlo ray tracing to accomplish RT. 
The simplicity of this approach allows new physics to be in- 
cluded in SPHRAY very easily and provides a framework that 
is conducive to parallelization. The accuracy of raytracing 
methods in general is determined by their ability to cover the 
simulation volume with a sufficient number of rays. This is a 
computationally expensive process and various approaches 
have been presented in the literature. 

Oxley & Woolfson (2003) combined a raytracing scheme 
with SPH in which the emitters and absorbers of radiation 
are at the deepest levels of a Barnes & Hut tree. In this 
scheme, the internal energies of the SPH particles are modi- 
fied by interpolating the changes in energy of the tree leaves 
onto the particles and vice versa. These radiative calcula- 
tions are made in between hydrodynamic time steps and 
allow for the coupling of the two processes. Stamatellos & 
Whitworth (2005) use a similar Barnes & Hut ray tracing 
scheme, but supplement the radiative transfer cells from the 
tree with a further star grid around sources. 

Both of these methods, in effect, propagate photon 
packets through a randomly chosen optical depth and then 



allow them to be absorbed or scatter. In this way, the tem- 
perature and emergent spectra can be calculated, however 
neither method solves for the ionization fraction explicitly 
but instead uses either a temperature-opacity or a total 
density-opacity relation. In addition, they raytrace through 
the adaptive grids generated from the Barnes & Hut trees 
and not the SPH particles themselves. 

Kessel-Deynet & Burkert (2000), taking advantage of 
the neighbor lists already in place in SPH simulations, in- 
troduced a fast method to find the optical depth from a 
source to a target particle. Variations of this method have 
since been utilized in Susa (2006) and Dale et al. (2007) 
to construct radiative transfer schemes in which each parti- 
cle influenced by the radiation of a source, becomes a tar- 
get of that source during the radiative update. A decision 
concerning how to treat the target particle is made based 
on the optical depth along the ray connecting the source 
and the target. Yoshida et al. (2007) have presented a re- 
lated method in which photon arrival times are calculated 
for each particle surrounding a given source by integrating 
the ionization front jump conditions along rays. After the 
photon arrival time for a source-particle pair, the photoion- 
ization rate is computed in the optically thin limit. Pawlik 
& Schaye (2008) have recently introduced an SPH / pho- 
ton packet based radiative transfer scheme in which packets 
are emitted into cones around the sources and propogated 
among the particles using the neighbor search lists. 

SPHRAY differs from the above methods in that it first 
uses a fast Axis Aligned Bounding Box (AABB) test to 
find the intersections of a ray and a Barnes & Hut tree 
that stores the particles. Next, the particles in these tree 
leaves are tested to see if they intersect the ray, yielding 
their impact parameter. In this way, every particle whose 
smoothing volume is intersected by the ray is stored in a 
ray list. In the course of moving along the ray from the 
source, the ionization and temperature state of each par- 
ticle is updated leading to more particle updates per ray. 
In this sense, SPHRAY shares the same ray- update ideas as 
the Monte Carlo ray tracing code CRASH by Maselli et al. 
(2003), but is applied to SPH density fields. With a sufficient 
number of rays, the native SPH resolution can be preserved. 

The format of this paper is as follows. In §2 we review 
the basic equations governing radiative transport and the 
evolution of the ionization and temperature state of a cos- 
mological gas. We also derive approximate analytic solutions 
necessary for an iterative numerical solution. In §3 we de- 
scribe the algorithms used by SPHRAY . In §4 we present the 
results of several standard RT tests. These tests were chosen 
to be the same as those in a recent radiative transfer compar- 
ison project (Iliev et al., 2006). They include: (1) isothermal 
HII region expansion; (2) HII region expansion with evolving 
temperature; (3) I-front trapping and shadowing by a dense 
clump; (4) multiple sources in a cosmological density field. 
We make our closest comparisons of SPHRAY with CRASH 
(the only other Monte Carlo code in the project), and a code 
described in Susa (2006) called RSPH (the only other SPH 
code in the project). Finally, in §5 we discuss future applica- 
tions and improvements of SPHRAY as well as summarizing 
and discussing our results. 
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2 BASIC PHYSICS - RADIATIVE TRANSFER, 
IONIZATION, AND TEMPERATURE 
EVOLUTION 

In this section we review the equations of 3D radiative 
transfer, the ionization and temperature equations for a 
cosmological gas, and the approximations that are made 
in SPHRAY . To facilitate comparison with other radiative 
transfer codes we review some of the other approximations 
which can be made. 



2.1 Notation 

In what follows, we use Roman numerals to indicate the 
ionization state of an element (H,He) in the standard way. 
Elements without Roman numerals refer to the nuclei of 
atoms (or all ionization states). A subscripted n refers to 
the number density of an element (or a specific ionization 
state of an element). A subscripted x refers to the ratio of the 
number density of a specific ionization state to the number 
density of all nuclei of that element. A subscripted y refers 
to the ratio of the number density of the subscripted species 
to the number density of H nuclei, for example, 



Xnell 
2/HeIII 



n-Hcii 

^HoIII 

n H 



(1) 
(2) 



2.2 Radiative Transfer Equation 

The 3-D radiative transfer equation in a frame comoving 
with the expansion of the Universe can be written (e.g. Nor- 
man et al., 1998), 



l dl v 
c dt 



ft- VI„ 



— (v— 

c \ dv 



(3) 



where e„ and the emission and extinction co- 

efficients respectively, H = a /a is the Hubble parameter, 
a = a/a e is the scale factor at time t divided by the scale 
factor at time t e (when the photons in the ray were emitted), 
and I v = J(x, n, v, t) is the specific intensity. 

For photons with a mean free path A m f p much less than 
the Horizon size c/H, the classical radiative transfer equa- 
tion is a valid approximation. 



c dt 



+ n • V/„ = t v — k v I v 



(4) 



This local approximation holds fairly well before the 
percolation stage of reionization when the growing ioniza- 
tion bubbles are still insulated from each other by the opti- 
cally thick IGM. Care must be taken once the majority of 
the IGM is reionized and becomes optically thin allowing 
photons to travel distances greater than the simulation box 
length. The effect of these background fluxes from outside 
the simulation volume must be taken into account, especially 
for high energy photons which have longer mean free paths 
and the potential to ionize and heat the IGM after being 
redshifted. The treatment of these non-local fluxes should 
be tailored to the specific problem at hand and so were not 
'hard-wired' into SPHRAY . For the test cases presented in 
§4 they were not necessary. 



Another caveat to using the classical equation, as ex- 
plained in Abel et al. (1999), is that it is only valid when 
\udl v /dv\ < I v and hence only for continuum radiation. 
However, the classical equation can still be used for line ra- 
diation if the redshifted absorption (photo-ionization) cross- 
sections are used when determining k„ . 

If e„ and k„ can be approximated as constant, a time 
independent RT equation can be used. 



n • V/„ — e v — K v Iv 



(5) 



This is a good approximation for individual SPH par- 
ticles over a sufficiently short time, however (as is also dis- 
cussed in Abel et al., 1999) it breaks down close to sources 
and allows the possibility of ionization fronts that travel 
faster than the speed of light. This can be quantified by 
examining the ionization front jump condition for a single 
point source ionizing a uniform density, constant tempera- 
ture, Hydrogen gas, 



dri 
dt 



~ ° H Jo 



n c nnxmidr 



(6) 



where, n is the distance to the ionization front from 
the source, N is the number of photons per second emitted 
by the source, and an is the recombination rate. An upper 
limit on the radius, r c within which the ionization front has 
a speed greater than c is, 



r c < 



N 



'iimuc 



(7) 



Within this region, use of the time independent equa- 
tion breaks down. In a raytracing scheme, this can be 
avoided by stopping rays once they have reached a distance 
d = ct on where t on is the amount of time the source has 
been on. The photons that were in the ray can be saved 
and traced from the stopping point once enough time has 
elapsed. In practice this is not always necessary. For exam- 
ple, the first test presented in §4 has r c /r s = 6.9 x 1CP 3 
where the Stromgren radius, r s = 5.4 kpc, 

In SPHRAY , the diffuse component of the radiation field 
is modeled using the on-the-spot (OTS) approximation, or 
as a set of many point sources and so for all calculations 
we can set e„ = along the ray, further simplifying the RT 
equation, 



dlv 
dr 



which has the analytic solution, 
Iu(r) = I„(r )e- T(r) 
where 



n x (r) a(v) dr 



k v dr 



(8) 



(9) 



(10) 



In principle, k v should include contributions from every 
process that removes photons from the ray under considera- 
tion (photo absorption, Thomson scattering, dust, etc.). For 
the tests presented here, we consider only photo absorption, 
however it would be straightforward to add terms to account 
for other processes. 
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2.3 Ionization Equations 

In this section we review the equations that determine the 
time development of the ionization fractions. They represent 
the contributions from photo-ionization, collisional ioniza- 
tion and recombination. Analytic and time averaged solu- 
tions in the case of constant rates are derived for use in an 
iterative solution scheme which relaxes the stringent con- 
straints on the time step. 

2. 3. 1 Chemistry 

SPHRAY follows the non-equilibrium evolution of six species 
[xHi,XHii,XHei,XHeii,XHein,i/e], only three of which are in- 
dependent. 

xm + xhii = 1 (11) 

XHcI + XHell + I'HcIII = 1 (12) 
2/e = 2/HII + 2/HeII + 2j/ H cIII + yz (13) 

where yz represents a constant background of free elec- 
trons from ionized metals. This background has a negligible 
effect on the evolution of the ionization fractions of H and 
He, but provides stability in the case of very small levels 
of ionization. We note that this is a subset of all atomic 
species relevant to primordial chemistry. Although inclusion 
of species involved in the formation of molecular Hydrogen 
[ H~, H2, Hj] is important in studies of primordial star for- 
mation, they have only a small impact on the evolution of 
the IGM. Primordial gas chemistry is discussed in detail in 
Anninos et al. (1997) and Abel et al. (1997). 



where, 

(— GhoI RhcII \ 

GhoI -(GHell + iW) ifaelll (22) 

GhoII — ^Helll / 

In general, every species with Ni B ionization states leads 
to an N- 1B x Ni s tridiagonal matrix. 

2.3.3 Analytic Solutions 

In order to proceed with a straightforward numerical inte- 
gration, the equations in the previous section are sufficient. 
However, the time steps are restricted by the stiff nature of 
the differential equations and so SPHRAY can also be run 
with an iterative ionization solver. This solver is based on 
the method used in the code C 2 -Ray presented by Mellema 
et al. (2006) where the detailed Hydrogen solution is given. 
For completeness we give the Helium solutions as well. The 
specific implementation in SPHRAY is described in §3. It re- 
quires time averaged analytic solutions for the ionization 
fractions which are derived below. 

If we assume that the Ga and Ri are constant the an- 
alytic solutions have the following form, 

xm(t) = x% + C],e vt (23) 
XHii(t) = x e ^ u + Cle vt (24) 



2.3.2 Differential Equations 

The processes that we will consider in the evolution of the 
ionization fractions are, recombination ojj, collisional ioniza- 
tion "/a, and photo-ionization Fa where A £ { HI, Hel, Hell 
} is one of the photo absorbing species and I G { HII, Hell, 
Helll } is one of the photo-ionized species. The equations 
can be written down directly, 



XHei(i) = x^ cl + Cu c e ll + C He e 



dxm 

dt 
dxmi 
dt 



dt 
dxacu 
dt 

dXHelU 

dt 



-GhiXhi + Rmixmi (14) 

GhiZhi — RauXnu (15) 

— GhcI^'HcI + -RHcIlZHell (16) 

GhcI^HcI — (GhcII + i?Hcii)a;Hoii + 

-RhcIII^HcIII (17) 

GhcII^'HcII — RhcUIXHcUI (18) 



where we have grouped the ionizing terms ( writing 
Ga = Ta +7a«. ) together, and included the electron num- 
ber density in the recombination term (Ri — a/n c ). In ma- 
trix form, 



x H 

XHo 



= M H x H 

= MhcXHc 



(19) 
(20) 



XHell(t) = xZ U + CLe Xlt +C^y 

XHain(i) = xZ m + Cke Xlt + C'Le X2t 
with the equilibrium solutions given by, 

eo Ran 



Ghi + Rmi 

Ghi 
Ghi + Rmi 



-RhoII-RhoIII 



ifecII-RHcIII + -RhcIIiGhoI + GhoiGhoII 

RhcUiGhoI 
RhcIiRhcUI + RhcUiGhcI + GHolGHell 

GhoiGhoII 
RhcIiRhcUI + RhcUiGhcI + GhoiGhoII 



(25) 
(26) 
(27) 



(28) 
(29) 

(30) 
(31) 
(32) 



Each system has one eigenvalue equal to zero corre- 
sponding to the equilibrium solutions. The non-zero eigen- 
values and the other constants can be expressed in terms of 
the Ga and Ri. For Hydrogen, 



G H 
G h 



Azhi = x m — X H i 
Axhii = a; H ii — %n 
— (Ghi + Rmi) 



(33) 
(34) 
(35) 
(36) 
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where a; HI and 

■^Hii *^r6 initial values and v is the non- 
zero eigenvalue for the Hydrogen system. For Helium, the 
expressions contain more terms, but they can be easily writ- 
ten down with the use of some notation. The matrix SHe of 
eigenvectors that will diagonalize Mh c is, 



(XHelll) = ^Helll" 



AiAt 



+ 



CL (e X2At - 1) 
\ 2 At 



(51) 



SHe = S'JL — 




(37) 



where Ai and A2 are the non-zero eigenvalues of the 
Helium system. 



Ai = -(s+p) 
A 2 = -(s-p) 

with 



= g (^HcII + RhcUI + GhoI + GHoIl) 

= RhcIiRhcUI + -RhoIIiGhcI + GhciGhcII 



d 



(38) 
(39) 



(40) 
(41) 
(42) 



We require the six constants G Ho in terms of the values 
Ga and Ri . Because the Helium system is constrained by the 
three differential equations and the fact that the ionization 
fractions must sum to one, there is some choice in the way we 
do this. SPHRAY follows XHeii and XHein and so a convenient 
relation is, 



GHe 

r 2 



= 6S 
= cS 

with, 



1.2 

He 
1,3 
lie 



Ghc — c 



G Ho — ^Hc 
Gwn = cSti' 



Azhciii — AxHeiiS H ' e 



3,2 _ g3,3 



3 Hc 

c 3,2 
I b He 



He 
AxhoIII 



-,3,2 
'He 



(43) 
(44) 



(45) 
(46) 



At the heart of this iterative method is the use of time 
averaged ionization fractions, optical depths, and photo- 
ionization rates to take larger time steps than would nor- 
mally be possible. The time averaged ionization fractions 
(x) = J Q x(t)dt of the above solutions are, 



2.4 Temperature Equation 

There are three terms in the temperature evolution equa- 
tion. One for the photo heating H, one for various atomic 
cooling 5 processes A, and one for the change in temperature 
due to the change in the number of free particles. 

dT 2 . s Tdn 

~dt 



"inks 



(n-A) 



n dt 



(52) 



2.4-.1 Photo Heating Term 

The term TL accounts for the kinetic energy of the photoion- 
ized electrons which quickly gets transferred to the other 
particle species. Let us suppose that a unit volume of gas 
absorbs 7V 7 photons per second. The fraction of absorption 
due to a given particle species is proportional to the optical 
depth of that species through the volume. The photo-heating 
rate for a monochromatic ray, H — Hm + HhcI + HucU can 
be simply expressed in this way, 

thi , 



Hm = N. 



7 (hv — hum) 

Tall 



HhoII = N 1 



N^—SL^hv — /lZ/Hel) 



Tall 
THell 
Tall 



(hl> - /lI/Hell) 



(53) 
(54) 
(55) 



where vm,vn e i, and ^Hcii are the ionization threshold 
frequencies. 

2.4.2 Atomic Cooling Term 

The atomic cooling function, A, includes the following phys- 
ical processes, 

• Collisional Ionization Cooling ( £ ) 

• Collisional Excitation Cooling ( tp ) 

• Recombination Cooling ( rj ) 

• Bremsstrahlung Cooling ( (3 ) 

• Compton Heating/Cooling ( \ ) 

The numerical value of A = A(ua, 
calculated using the rates detailed in Appendix A. 



(xm) = x HI + — {e - 1, M 

I \ eg . G H / „Ai , \ 1 



At' 



(XHel) = ^H g el + 



Gfjc (e AlAt - 1) G He (e^ A< - l) 



XiAt 
(XHell) = x% eII + 



\ 2 At 



G Hc (e^ - 1) G Hc (e^ A< 1) 



XiAt 



+ ■ 



A 2 At 



(47) 
(48) 

(49) 
(50) 



3 NUMERICAL TECHNIQUES 

In this section we outline the numerical techniques used to 
solve for the ionization and temperature state. SPHRAY is a 
Monte Carlo code and is based on sampling the radiation 
field along 1-d characteristics. This is accomplished by trac- 
ing rays that extend a predefined length. The length can be 
chosen in a number of ways. For vacuum boundary condi- 
tions the obvious choice is to terminate the rays at the edge 

5 note that Compton scattering can have a negative contribu- 
tion to the cooling function if the temperature of the background 
radiation field is greater than the gas kinetic temperature 
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of the volume. For reflective or periodic boundary condi- 
tions a criterion can be applied to the properties of photons 
in the ray (for example when the flux has dropped below a 
specified value) or a hard limit on the length of a ray can 
be set (for example 2 box lengths). The impact parameter 
for every particle-ray intersection is calculated using a fast 
AABB test allowing for the calculation of the photon flux 
at all the particle-ray intersections. 

3.1 Sources 

Any point in the simulation volume can be specified as the 
beginning of a ray. Currently, SPHRAY is configured to treat 
point sources whose properties are specified by an input file, 
recombination rays from ionized SPH particles, and back- 
ground fluxes by specifying points on the simulation volume 
walls as sources. 

3.1.1 Diffuse Recombination Radiation 

Here we review the recombination processes that produce 
ionizing photons in a H/He gas (e.g., Osterbrock, 1989). The 
following free-bound transitions produce continuous spectra. 



HII + e^HI(l 2 S)+7 («13.6eV) (56) 

HeII + e^ HcI(l 1 S)+7 (^24.6eV) (57) 

Helll + e -> HeIl(2 2 S) + 7 (wl3.6eV) (58) 

Helll + e -> HelI(2 2 P) + 7 (wl3.6eV) (59) 

Helil + e -> HeII(l 2 S) +7 (^54.4eV) (60) 



These spectra can be calculated exactly using the Milne 
relations. For Hydrogen, the emission coefficient for the 
above process is (Osterbrock, 1989), 

£H = — ) **< " (61) 

where is the photo absorption coefficient of Hydrogen in 
the ground state. Similar relations can be derived for the 
other free-bound processes, however in practice these highly 
peaked spectra can be approximated by delta functions just 
above the appropriate threshold (in parentheses in the 
equations above). 

Following free-bound captures to excited Helium states, 
the following bound-bound transitions can also produce ion- 



izing photons. 

HeI(2 3 S) -> Hel^S) + 7 19.8 eV (62) 

HeI(2 1 P) -> Hel^S) + 7 21.2 eV (63) 

HeI(2 1 S) -> Hel^S) + 2 7 E20.6 eV (64) 

HeII(2 2 P) -> HeII(l 2 S) + 7 40.8 eV (65) 

HeII(2 2 S) -> HeII(l 2 S) + 2 7 E40.8 eV (66) 



The various bound-bound transitions above have relative 
probabilities that depend on the environment (free electron 
density, temperature, ionization state) and so the weights to 

6 It is a coincidence that electron captures by HcIII directly to 
the n=2 level of Hell have a spectrum peaked at the Hydrogen 
threshold. 



give to these processes should be tailored to specific appli- 
cations. 

The most straight forward way to deal with this diffuse 
radiation is to use the on-the-spot (OTS) approximation. 
The OTS approximation makes the assumption that recom- 
bination photons are absorbed in the vicinity (the same SPH 
particle) of the point where they are emitted. Computation- 
ally, this means no ray tracing is necessary for these pho- 
tons. For pure Hydrogen simulations, the OTS approxima- 
tion amounts to using the reduced recombination rates in 
the appendix (case B) 7 . Here, one is making the assump- 
tion that each electron capture directly to the ground state 
produces a photon that ionizes a nearby Hydrogen atom and 
the two actions effectively cancel one another. 

For simulations involving Helium, using the case B rates 
would be making the assumption that each Hell recombina- 
tion to the ground state ionizes a nearby Hel atom while each 
Helll recombination to the ground state ionizes a nearby 
Hell atom. This is a simple first approximation, but some 
of the Hell ground state captures will ionize HI, while some 
of the Helll ground state captures will ionize HI and Hel. 
To account for this would require a more detailed adjust- 
ment of the recombination and photoionization rates for the 
species involved. The most computationally intensive option 
is to trace rays for each of these recombination processes and 
thereby take account of the fact that some of the photons 
will not be absorbed in the SPH particle where they were 
created. Again, the level of detail used should be guided 
by the application at hand. With SPHRAY it is possible to 
choose either the recombination ray or the OTS approach. 
We plan to explore the accuracy of the OTS approximation 
in various geometries and densities in future work. 

3.2 Optical Depth in SPH 

Ray tracing solutions to the radiative transfer problem solve 
the equation along ID characteristics. As such, an estimate 
of the optical depth along these characteristics is central 
to the problem. In the SPH formalism, a continuous den- 
sity field is represented by a number of discrete fluid ele- 
ments (particles) with smoothing lengths hi. These smooth- 
ing lengths 8 are usually defined to keep a constant mass 
M S ph inside the smoothing volume Vi = |7i7if. The prop- 
erties of the fluid at any point are then estimated by aver- 
aging over all N particles in the simulation weighted by a 
smoothing kernel. In practice, one only averages over nearby 
particles, but this definition is equally valid and useful in the 
derivation to follow. As an example, the density p(r\) at the 
position 17 of the i th particle is estimated as, 

JV N 

p(n) K^mjWQn - r 3 \,h) =^m j W(r ij ,h) (67) 

where rrij is the mass of the j th nearest particle and 
W(rij,h) is a smoothing kernel. 



7 Case A rates refer to recombinations to all atomic levels. Case 
B rates refer to recombinations to all but the first atomic level so 
that ajj = aj| — . 

8 Throughout this work we use the convention that the smooth- 
ing kernel goes to zero at h and not 2h 
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An estimate of the fluid property need not be made at 
the position of a particle. An averaged value for any fluid 
property can be defined for an arbitrary point in space using 
two techniques. One is the "scatter" method in which the 
desired quantity is calculated by averaging over every parti- 
cle whose smoothing volume includes the point in question. 



P(ri) 



(68) 



For this case, a different smoothing length hj is used in the 
kernel W for each term. In the "gather" method, a smooth- 
ing length is defined for the arbitrary point and the desired 
quantity is calculated as an average over the particles within 
this smoothing length. 



JV 

p( r ~^2mjW(r i3 ,hi) 

j'=i 



(69) 



For this case, only one smoothing length hi is involved. For 
definiteness, a popular spline kernel (Monaghan and Lat- 
tanzio 1985) is, 



W(r,h) 



8 
ivh 3 



1 -6(tt) 2 + 6(^) 3 0<r<| 

2(1 -^f f i<r<h (70) 

r > h 



In our ray tracing scheme, we would like an estimate 
of the column depth N c d along a ray that has intersected a 
number of SPH particles. Formally this is, 



iVcd 



/ p(r)d/« / ^m,W/(n 3 A)dZ (71) 
Jo Jo j=1 



where the scatter interpretation has been used. The sum- 
mation extends over all particles and I parameterizes the 
distance along the ray. Interchanging the order of integra- 
tion and summation we have, 



N, 



cd 



i=l J ° 



rrijW(rij, hj)dl 



(72) 



This amounts to a line integral through the smoothing 
kernel of each particle whose smoothing volume is pierced 
by the ray. This integral is calculated by tabulating it as a 
function of impact parameter b for one value of h, namely 
h = 1. We can recover the line integral for any b and h 
through a rescaling of the tabulated value. This technique 
delivers the optical depth to a point along the ray, as well 
as the contribution from each particle in the raylist to that 
optical depth. 

There are two related approximations that come into 
play here. The first involves the reordering of the terms so 
that those involving the same particle are adjacent. This 
is valid as long as the density doesn't vary much within a 
smoothing length which is true by construction. The next 
involves the endpoint of the ray. Given this reordering, the 
contribution to the density (at the terminus of the ray) 
from particles which have centers further down the ray but 
smoothing volumes' that contain the end point, will be un- 
accounted for. This is a small correction for the reason given 
above, and the fact that the photons are usually nearly all 
absorbed when the ray is terminated. 



3.3 Photon Packets 

In our Monte Carlo method, the radiation field is discretized 
into photon packets and transported along rays. Each packet 
contains a large number of monochromatic photons sam- 
pled from an arbitrary spectral energy distribution (SED). 
The direction along which a photon packet is transported 
is also sampled from an emission profile distribution. This 
is true whether the packet represents emission from a point 
source, diffuse recombination emission, or background ra- 
diation originating outside the computational volume. The 
starting locations for rays can be any point within the com- 
putational volume including points on the faces of the sim- 
ulation volume. 

For each ray that is cast, a source is selected at ran- 
dom weighted by its luminosity. This ensures a population 
of photon packets with roughly the same energy as opposed 
to the same number of photon packets being traced from 
each source regardless of their luminosities. 

The base resolution of a simulation is determined by 
how many SPH particles N p are used to sample the contin- 
uous density field. The degree to which the sampling of the 
radiation field approaches the base resolution is determined 
by the number of rays traced N T . For isotropic sources and 
homogeneous density fields, the average number of particle 
intersections per ray is ~ A^ 3 . It follows that the total 

1/3 

number of intersections is on the order of N T N P ' and that 
the average number of intersections per particle is N T N P 2 ^ 3 . 
This gives a rough estimate of how many times the radiation 
field is sampled at each particle. In practice, particles closer 
to sources will be sampled more often and it is better to 
conduct a convergence study than to rely on pre-calculated 
estimates of resolution. 



3.4 Particle-Ray Intersections 

Once a photon packet has been constructed, it is propagated 
along a ray. SPHRAY uses a data object called a raylist to 
store the intersections of a ray drawn from the source, and all 
the SPH particles whose smoothing volumes are pierced by 
the ray. This can be done using vacuum, or periodic bound- 
ary conditions (in the latter case a maximum length must 
be specified). The search for these intersections needs to be 
as efficient as possible. 

Because the smoothing lengths within a cosmological 
box can vary by more then three orders of magnitude, we or- 
ganize the particles into an oct-tree. This manner of storing 
particles is common and many SPH codes produce an oct- 
tree during the course of hydrodynamic calculations. Our 
code could be trivially modified to use a pre-constructed 
tree although currently SPHRAY constructs its own for each 
density field it ray traces. We augment the standard oct-tree 
by associating with every cell an axis-aligned bounding box 
(AABB) that is just large enough to encompass the smooth- 
ing lengths of all the particles in that cell (Figure 1). 

Each cell of our tree contains either particles or daugh- 
ter cells. The maximum number of particles in a leaf is spec- 
ified in a configuration file. We have found 12 to be a rea- 
sonable choice. The search for particle-ray intersections pro- 
ceeds exactly as in the case of a simple SPH neighbor search. 
Starting with the root cell, the AABB of the cells are tested 
for intersection with the ray. In case of intersection, the cell 
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Figure 1. Cartoon 2D representation of several SPH particles. 
The left panel shows the positions of the axis aligned bounding 
boxes and the right panel indicates the quad-tree cell boundaries. 



is opened and the search proceeds on the daughter cells. 
This process continues until a leaf with no further refine- 
ments is encountered, in which case the particles in the leaf 
are tested to see which of them are intersected by the ray. 
This also produces the impact parameters of each particle 
in the ray list. 

The intersection test of the ray with the AABB is done 
using pliicker coordinates (Mahovsky & Wyvill, 2004), a 
very fast method used in computer graphics. This method 
tests the ray against the edges comprising the silhouette of 
the AABB instead of testing against the individual faces, 
it is division free and consists of a number of simple dot- 
product operations. 



3.5 Solvers 

Once we have a photon packet and a list of the ray-particle 
intersections stored in a raylist, we can proceed to update 
all the particles in the raylist. SPHRAY offers two choices 
for this task. The first is an adaptive Runge-Kutta (RK) 
method. A set of formulas due to Fehlberg (1970) provide 
solutions that are accurate to fifth order in the time step. 
Step size control is provided using the truncation error as 
estimated by the embedded fourth order solution. The Cash- 
Karp coefficients Cash & Karp (1990) are used to take the 
variable length time steps. This option is included because it 
is simple and could be easily modified if a user needs to add 
extra physics into the solution routine. It is also guaranteed 
to conserve photons to an arbitrary accuracy by forcing the 
number of ionizations in a particle to equal the number of 
photons removed from a packet. 

The second solution method is an iterative solver based 
on time averaged photoionization rates and optical depths. 
The time averaging removes the need for very short time 
steps, but requires approximate analytic solutions and in 
the case of extremely long time steps will not conserve pho- 
tons exactly. The main advantage is that the number of it- 
erations necessary to obtain a converged solution can be 
much smaller than the number of RK time steps necessary 
to obtain the same solution. This method was introduced in 
Mellema et al. (2006) and a more detailed description of it 
can be found there. 



3. 5. 1 Runge-Kutta 

For each intersection we determine the time tn since the 
particle has last been intersected by a ray. The photon flux 
N y at the particle is estimated as the photons left in the ray 
Ni divided by tn. This is taken to be the first guess for a 
time step in the solution of the system of coupled differential 
equations (Eqs. 15,17,18, and 52). The optical depth through 
the particle is estimated using the technique described in 
§3.2. This allows the calculation of the total photoionization 
rate T which along with the values of Ga and Ri are all that 
is necessary to calculate the right hand sides of the equations 
mentioned above. The photoionization rate is 



7BH 

Mr, 



Y 

Xxm + -j( XHeI + ^Hell, 



(73) 



where JV 7 is the photon flux at the particle, At is the 
optical depth through the particle, mu is the mass of a Hy- 
drogen atom, M p is the mass of the particle, X is the Hydro- 
gen mass fraction and Y is the Helium mass fraction. Here, 
xh b i and :EHeii should be set to zero for frequencies less than 
their respective thresholds. T for the individual species is 
calculated from the ratio of their optical depths to the total 
optical depth. 



At 



(74) 



The number of photons absorbed by a particle iV a is 
obtained by solving an extra differential equation, as the 
photoionization rate is allowed to vary for each substep that 
the RK routine takes. 

The photon packet is followed along a ray until the frac- 
tion of photons left is below a threshold or until the photon 
packet has reached a pre-determined distance along the ray. 
A typical value for the photon tolerance is 1.0 x 10~ 10 of 
the initial photons in the packet. This determines the level 
of photon conservation and can be set arbitrarily low. 



3.5.2 Iterative 

If the iterative solution method is chosen, the default time 
step tn can be used for most updates. The first step in this 
solution method is to initialize the time averaged ionization 
fractions to the current values in a particle. These are used 
to make a first guess at the time averaged optical depth and 
photoionization rate. 



<r> = at t ( 



M v 



l-e-< A ^>) X 
Y 

X (xm) + — ((sHel) + (zHdl)) 



(75) 



This time averaged photoionization rate is used to cal- 
culate the time averaged ionization fractions (Eqs. 52-56) 
which are in turn used to update the time averaged opti- 
cal depths. The optical depths can then be used to find a 
new photoionization rate and the iteration proceeds until 
we have reached convergence in the electron number density 
and temperature. Because the heating and cooling rates are 
themselves functions of temperature (as opposed to the re- 
combination and collisional ionization rates which are not 
functions of the ionization fraction), the temperature is al- 
ways updated using the RK routine and the ionization state 
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at each iteration. Once the iterations have converged, the 
final state of the particle is calculated using eqs. 28 - 32. 



4 CODE VERIFICATION 

Here we present the results of SPHRAY on the tests outlined 
in the radiative transfer comparison project paper by Iliev 
ct al. (2006). We will describe the tests briefly, but refer 
the reader to the reference for details. For tests 1 and 2, 
the initial conditions were set up using SPH particles that 
were evolved into a glass state using the code GADGET- 
2 (Springel, 2005). The distribution of particles extended 
two smoothing lengths further then the required box size in 
order to avoid edge effects in the density field and vacuum 
boundary conditions were used. 

Tests 1 and 2 contained exactly 128 3 SPH particles 
within the 6.6 3 kpc 3 volume while test 3 contained 2,135,842 
particles in this same volume. For the third test, 853,442 
particles were first evolved into a glass and then SPH parti- 
cles were randomly placed inside the clump until the density 
there had reached 200 times the density outside the clump. 
The initial conditions for the cosmological test are discussed 
in §4.1.4. The SPH density field as well as the temperature 
and ionization fraction variables have been interpolated onto 
a 128 3 grid to make the surface plots and for submission to 
the Comparison Project website and so for those figures we 
will refer to grid cells. 

4.1 Test 1. Pure hydrogen isothermal H II region 
expansion 

The first test considers the growth of a Stromgren sphere 
in a uniform density field consisting of pure hydrogen. 
The source, placed in the corner of a 6.6 kpc box, emits 
iV 7 = 5.0 x 10 48 Rydberg photons per second. The density 
of hydrogen is nn = 1.0 x 10~ 3 cm~ 3 , and the temperature 
is fixed at T = 10, 000 K. The system is evolved for 500 Myr 
or approximately four recombination times. The state of the 
system is examined at 10 and 100 Myr when the I-front is 
growing quickly and at 500 Myr when the photoionizations 
and recombinations have balanced and the HII region has 
reached its final Stromgren radius. 

This simple test has the advantage that numerical re- 
sults can be compared directly with an analytic solution 
for the position and velocity of the I-front versus time. 
SPHRAY finds agreement with these solutions at the several 
percent level (Figure 3). The analytic front width, defined 
as the distance over which the neutral fraction goes from 
xm = 0.1 to xm = 0.9 is nt ~ 18A m f p « 14 grid cells (Iliev 
et al., 2006) and is also reproduced by SPHRAY . This can be 
seen in figure 2, by noting that the contours are at xhi = 0.1 
and xm = 0.9. 

The size of the highly ionized proximity region produced 
by all the codes in the comparison project can be seen by ex- 
amining the neutral fraction lines in figure 3. Our solution 
follows closely that produced by RSPH. There are several 
codes that produce slightly smaller proximity regions then 
the rest. These codes may not have converged or be un- 
able to reach the native resolution of the density field as 
SPHRAY produced similar results before converging as can 
be seen in figure 4. 



The anisotropies seen in the surface plots of SPHRAY 's 
results have three major sources: the Monte Carlo method 
used, Poisson fluctuations in the sampling of solid angle by 
the rays, and fluctuations in the density field due to the 
SPH glass. In future work, it may be useful to use a low dis- 
crepancy sequence instead of uniformly distributed pseudo 
random numbers to generate ray directions, however given 
the fact that these anisotropies can be reduced by tracing 
more rays and the agreement of the radially averaged profiles 
we consider this a minor problem. 

4.2 Test 2. HII region expansion: the temperature 
state 

Test 2 is identical to Test 1 except the gas temperature 
is now allowed to vary, starting from an initial value of 
T = 100K and the spectrum of the ionizing source is taken 
to be that of a T = 10 5 K blackbody. This is a more de- 
manding test as the heating times are very short in compar- 
ison with the other characteristic times and now the rays 
must sample frequency space as well. Multi-frequency pho- 
ton packets would alleviate the need for this extra sampling 
and are planned as an additional option in future versions 
of SPHRAY . 

The radially averaged ionization ( figure 7 ) and tem- 
perature ( figure 8 ) profiles that SPHRAY produced for 
this test again follow very closely those produced by RSPH. 
SPHRAY correctly produces a preheated region ahead of the 
ionization front due to the accurate treatment of high en- 
ergy photons with long mean free paths and a large amount 
of energy to deposit as heat. 

4.3 Test 3. I-front trapping in a dense clump and 
the formation of a shadow 

In this test, a cold dense spherical clump of hydrogen gas is 
embedded in a hot diffuse background. The dimensions of 
the simulation box are the same as above (6.6 3 kpc 3 ). The 
clump has a radius r c = 0.8 kpc and its center is located at, 
x c = (5.0, 3.3, 3.3) kpc. The density contrast is ni n /n out = 
200 with n ou t = 2 x 10~ 4 cm~ 3 . The gas in the clump is 
set to a temperature of T; n = 40 K while the gas outside 
the clump is initialized to T out = 8000 K. The entire x = 
side of the simulation box is taken to be a T = 10, 000K 
blackbody source with a constant photon flux, F = 10 6 s _1 
cm~ into the box. The test is designed so that the initially 
fast moving ionization front will be trapped in the clump 
due to the higher recombination rate there. 

In figures 9 and 10 we show the ionization and tem- 
perature profiles along a small cylinder through the axis of 
symmetry. For the Comparison Project grid data, we used 
the four central columns of grid cells (where each grid cell 
is « 0.05 kpc in length) and for our SPH data we used all 
the particles whose centers lie within 0.05 kpc of the axis of 
symmetry. The variation between all the codes is larger for 
this test than the first two. SPHRAY produces results that 
lie between those of RSPH and CRASH for the ionization 
profiles and results that follow RSPH for the temperature 
profiles. 

This test features many of the effects that contribute 
to a reprocessing of the intergalactic background radiation 
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Figure 2. Comparison Project Test 1 (HII region expansion in a gas at constant density and temperature). Surface plot of xjji cut 
through the simulation volume at coordinate z = at t = 10 (left), 30 (middle), and 500 (right) Myr. The contours are at xjji =0.1 and 
^HI = 0-9- The number of rays traced N t by SPHRAYis q X 10 s where q = f/500 is the fraction of elapsed time. The axes are measured 
in grid cells. 




Figure 3. Comparison Project Test 1 (HII region expansion in a gas at constant density and temperature), xhi and xhii radial profiles 
at t = 10 (left), 100 (middle), and 500 (right) Myr. The number of rays traced N t by SPHRAYis q X 10 7 where q = i/500 is the fraction 
of elapsed time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple line), RSPH (dashed green line), 
C 2 -Ray (small dashed light blue line) and the results from all the other codes in the comparison project that completed this test (dotted 
orange lines). 



field. Specifically, self shielding, shadowing, and spectral 
hardening. SPHRAY produces an ionization front that moves 
quickly to the clump and is trapped near the center as it 
should be. The differences with the other codes mostly have 
to do with the amount of shadowing and self shielding (ion- 
ization and temperature). The high energy photons that are 
sampled by SPHRAY produce slightly weaker shadows di- 
rectly behind the clump than most of the other codes. The 
results are much more pronounced for the temperature than 



they are for ionization. In the last panel of figure 11 the self 
shielded section of the clump and the area in the shadow 
of the clump are only slightly ionized, however in the last 
panel of figure 12 the high energy photons have raised the 
temperature of the whole initially cool clump to w 10, 000 
K and begun to heat the shadowed region behind it above 
its initial temperature of 8, 000 K, 



Radiative Transfer with SPHRAY 11 




Figure 4. Comparison Project Test 1 (HII region expansion in a gas at constant density and temperature), xhi and xhii radial profiles 
at t = 10 Myr with N% = q X 10 6 (left), JVt = q X 10 7 (middle), and 7V t = q X 10 s (right) where Nt is the number of rays traced and 
q = 10/500 is the fraction of elapsed time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple line), 
RSPH (dashed green line), C 2 -Ray (small dashed light blue line) and the results from all the other codes in the comparison project that 
completed this test (dotted orange lines). 





Figure 5. Comparison Project Test 2 (HII region expansion in a gas at constant density with varying temperature). Surface plot of 
xm cut through the simulation volume at coordinate z = and t = 10 (left), 100 (middle), and 500 (right) Myr. The contours are at 
xhi = 0.1 and xm = 0.9. The number of rays traced JVt by SPHRAY is q X 10 s where q = t/500 is the fraction of elapsed time. The axes 
are measured in grid cells. 



4.4 Test 4. Multiple sources in a cosmological 
density field 



The most realistic test preformed in the RT comparison 
project involved sources placed in the 16 most massive halos 
of a cosmological simulation snapshot at z=9. The box size 
is 0.5 ft -1 comoving Mpc and the gas was initially neutral 
with a temperature of 100 K. The sources were assumed to 



emit / 7 = 250 photons per atom over t a = 3 Myr leading to 
a photon flux iV 7 of, 



JV 7 = /, 



(76) 



where M is the halo mass. The system was then evolved for 
0.4 Myr. 

In order for SPHRAY to complete this test it was neces- 
sary to convert the density field data from a grid based rep- 
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Figure 6. Comparison Project Test 2 (HII region expansion in a gas at constant density with varying temperature). Surface plot of the 
temperature cut through the simulation volume at coordinate z = and t = 10 (left), 100 (middle), and 500 (right) Myr. The number 
of rays traced 7V t by SPHRAYis q X 10 8 where q = i/500 is the fraction of elapsed time. The axes are measured in grid cells. 
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Figure 7. Comparison Project Test 2 (HII region expansion in a gas at constant density with varying temperature). Spherically averaged 
xhi and xjjii profiles at t = 10 (left), 100 (middle), and 500 (right) Myr. The number of rays traced At by SPHRAYis q X 10 7 where 
q = i/500 is the fraction of elapsed time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple line), 
RSPH (dashed green line), C 2 -Ray (small dashed light blue line) and the results from all the other codes in the comparison project that 
completed this test (dotted orange lines). 



resentation into an SPH density field. This was accomplished 
using the following procedure: we start from a smooth glass 
distribution representing a constant density field equal to 
the peak density p ma x of the grid data. Particles are then 
selected for removal according to a comparison of the den- 
sity p p at the particle position with a random density px 
where < px < Pmax- This gives a density field which is 
smoother in high density regions than it would be if a nor- 
mal acceptance-rejection test with random particle locations 



was used. In low density regions the noise is higher because 
the distribution becomes less glass like as particles are re- 
moved. The result was a close approximation to the gridded 
density field using 1,957,344 SPH particles. The test does 
show some extraneous noise from the conversion to parti- 
cles, but this noise is smaller than the difference between 
different methods so for our comparison, the simple proce- 
dure outlined above was good enough. SPHRAYis intended 




Figure 8. Comparison Project Test 2 (HII region expansion in a gas at constant density with varying temperature). . Spherically 
averaged temperature profiles at t = 10 (left), 100 (middle), and 500 (right) Myr. The number of rays traced Nt by SPHRAYis q X 10 7 
where q = t/500 is the fraction of elapsed time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple 
line), RSPH (dashed green line), C 2 -Ray (small dashed light blue line) and the results from all the other codes in the comparison project 
that completed this test (dotted orange lines). 




Figure 9. Comparison Project Test 3 (I-front trapping in a dense clump), xjji and xhii profiles along the axis of symmetry at t = 1 
(left), 3 (middle) , and 15 (right) Myr. The number of rays traced JVt by SPHRAYis q X 10 8 where q = t/15 is the fraction of elapsed 
time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple line), RSPH (dashed green line), C 2 -Ray 
(small dashed light blue line) and the results from all the other codes in the comparison project that completed this test (dotted orange 
lines). 



to be run on density fields produced by SPH hydrodynamic 
simulations in which this sort of noise would not occur. 

The results from SPHRAY are compared with those of 
C 2 -Ray (Mellema et al., 2006), FTTE (Razoumov & Cardall, 
2005), and CRASH (Maselli et al., 2003) in figures 13 to 16. 
General features of the ionization field including the extent 
and shape of the ionization front, the shadows from dense 
clumps, and the shape of the neutral island near the center 
of the slice, are similar in all codes. Although SPHRAY and 
CRASH share the most similarities, including the sampling 
of high energy photons, we see some slight differences in the 



peak ionization and temperature values they produce in the 
central part of the highly ionized region. 



5 DISCUSSION 

We have presented the cosmological SPH raytracing code 
SPHRAY and discussed the results of several radiative trans- 
fer problems by way of validation. SPHRAY employs a Monte 
Carlo approach to ray tracing applied directly to the SPH 
particle distribution native to a large fraction of current 
astrophysical and cosmological hydrodynamic simulations. 




Figure 10. Comparison Project Test 3 (I-front trapping in a dense clump). Temperature profiles along the axis of symmetry at t = 1 
(left), 3 (middle), and 15 (right) Myr. The number of rays traced JVt by SPHRAYis q X 10 s where q = t/15 is the fraction of elapsed 
time. Shown are the results from SPHRAY (solid red line), CRASH (dash-dot-dot-dot purple line), RSPH (dashed green line), C 2 -Ray 
(small dashed light blue line) and the results from all the other codes in the comparison project that completed this test (dotted orange 
lines) . 




Figure 11. Comparison Project Test 3 (I-front trapping in a dense clump). Surface cut of the neutral fraction at t = 1 (left), 3 (middle), 
and 15 (right) Myr. The number of rays traced Nt by SPHRAYis q X 10 s where q = t/15 is the fraction of elapsed time. Shown are the 
results from SPHRAY 



The column density sums at the heart of the radiative trans- 
fer calculation are carried out using the SPH kernels so that 
no regridding of the data is necessary, maintaining the adap- 
tive Lagrangian nature that makes SPH attractive in the 
first place. The statistical nature of the Monte Carlo method 
makes the inclusion of arbitrary source spectra and emission 
profiles very straightforward. The simplicity of its implemen- 
tation also allows the future addition of more complicated 
physics as well as parallelization. However, a large number 
of rays must be traced to get a fair representation of the 



underlying probability distribution functions being sampled 
and to maintain angular resolution. This translates to the 
numerical problem of finding the intersection of numerous 
rays and spheres as quickly as possible. In order to do this 
we have applied a variant of the neighbor search techniques 
using oct-trees and a fast box-ray intersection test adapted 
from computer graphics, resulting in an efficient adaptive 
ray tracing code applicable to current astronomical hydro 
simulations. 

There are, in general, no analytic solutions to the 
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Figure 12. Comparison Project Test 3 (I-front trapping in a dense clump). Surface cut of the temperature at t = 1 (left), 3 (middle), 
and 15 (right) Myr. The number of rays traced JVt by SPHRAY is q X 10 8 where q = t/15 is the fraction of elapsed time. Shown are the 
results from SPHRAY 
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Figure 13. Comparison Project Test 4 (Multiple sources in a 
cosmological density field). Surface cut of the neutral fraction 
through the middle of the simulation volume at t = 0.05 Myr. The 
number of rays traced JVt by SPHRAY is q X 10 8 where q = t/0.4 
is the fraction of elapsed time. Beginning in the top left corner 
and proceeding clockwise are the results from C 2 -Ray, FTTE, 
SPHRAY, and CRASH. 
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Figure 14. Comparison Project Test 4 (Multiple sources in a cos- 
mological density field). Surface cut of the temperature through 
the middle of the simulation volume at t = 0.05 Myr. The number 
of rays traced JVt by SPHRAY is q X 10 8 where q = t/OA is the 
fraction of elapsed time. Beginning in the top left corner and pro- 
ceeding clockwise are the results from C 2 -Ray, FTTE, SPHRAY , 
and CRASH. 
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Figure 15. Comparison Project Test 4 (Multiple sources in a 
cosmological density field). Surface cut of the neutral fraction 
through the middle of the simulation volume at t = 0.2 Myr. The 
number of rays traced JVt by SPHRAYis q X 10 s where q = t/0.4 
is the fraction of elapsed time. Beginning in the top left corner 
and proceeding clockwise are the results from C 2 -Ray, FTTE, 
SPHRAY, and CRASH. 
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Figure 16. Comparison Project Test 4 (Multiple sources in a cos- 
mological density field). Surface cut of the temperature through 
the middle of the simulation volume at t = 0.2 Myr. The number 
of rays traced N t by SPHRAYis q X 10 s where q = t/0.4 is the 
fraction of elapsed time. Beginning in the top left corner and pro- 
ceeding clockwise are the results from C 2 -Ray, FTTE, SPHRAY , 
and CRASH. 



types of radiative transfer problems that occur with sources 
embedded in 3D density fields. Therefore we validated 
SPHRAY using the tests chosen by the Radiative Trans- 
fer Comparison Project. There is good agreement of 
SPHRAY results with codes that treat the same level of 
physics. 

We present the source code for SPHRAY on a compan- 
ion website 9 , together with a users guide and some example 
input snapshot files. An early version of this RT approach 
was used to study the structure of neutral hydrogen in the 
Universe at the time of reionization in Croft & Altay (2007) . 
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APPENDIX A: PHOTO-IONIZATION 
CROSS-SECTIONS AND ATOMIC 
COOLING/HEATING RATES 

The rates below make use of the following notation, 



A. 



J T 



(Al) 



where the Ta are the ionization energies of the photo 
absorbing species in temperature units, 



Thi = 157, 809K 
T Hc i = 285, 335K 
Thoii = 631.515K 
and, 

T 



(A2) 
(A3) 
(A4) 

(A5) 



10* K 

• Recombination Rates - Case A (Hui & Gnedin, 1997) 
fern 3 s" 1 ! 



Ohh = 1-269 x 10~ 13 



A 



1.0 + 



„ A Q O 1n -14 .,0.654 

«Hoii = 3.0 x 10 A HcI 



OHem = 2.538 x 10 -13 



(A6) 

(A7) 
(A8) 



• Recombination Rates - Case B (Hui & Gnedin, 1997) 
[cm 3 s _1 l 



ai u = 2.753 x 10~ 14 



1 1.500 
^HI 



http: / /www. sphray.org 



IQ+ ( A HI )' 
±.u T V 2.740/ 



(A9) 
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"Hell = i- 26 X 10 A HcI 



afein = 5-506 x 10" 14 



\ 1.500 
A HcII 



1.0 + 



'A H cIl\ - 407 



• Collisional Ionization Rates (Cen, 1992) [cm 3 s 1 ] 

5.85 x lO" 11 ^ _ Tm/T 
THi = = e m/ 

i + Vn 

2.38 x lO" 11 ^ _T HeI /T 
7Hei = = e HtI/ 

i + Vn 

= 5.68 x 10- 12 V^ c _ THaIl/T 

i + Vt¥ 

• Collisional Ionization Cooling (Cen, 1992) 
[ergs cm - ' 1 s" 1 ] 



(A10) 
(All) 

] 

(A12) 
(A13) 
(A14) 



. 1.27 x lO" 21 ^ _ Tm/T 

. 9.38 x 10" 22 ^ _T HeI /T 
CHei = ; — ; — t= e HeI/ n c n Ho i 



, 4.95 x 1Q- 22 V% -t Ho „/t 

ChcII = = e Hell/ n n He ii 

1 + V 

• Collisional Excitation Cooling 
1992) [ergs cm" 3 s" 1 ] 

7.5 X 10 19 -118348/To 
1 + y -L 5 
, 9.10 X l -27 r -0.1687 _ 

V'Hei = = e °n e n H eii 

1 + vTr, 

5.54 X 10 17 T 0,397 _473638/Tn 

1 + VT5 

• Recombination Cooling - Case A (Hui & Gnedin, 1997) 

[ergs cm" 3 s" 1 ] 



(A15) 

(A16) 

(A17) 
(Cen, 

(A18) 
(A19) 
(A20) 



?7hii = 1-778 x 10" 29 



rp \ 1.965 



i-o + (i^) <: 



n e nmi (A21) 



VHell = fcfcTbaHell^eriHell 



r&an = 1.4224X10" 28 



To A j 



i-o+ (M) c 



• Recombination Cooling - Case B (Hui & Gnedin, 1997) 



(A22) 

1.923 n e^HeIII ( A23) 



[ergs cm 3 s 1 ] 



rp ,1.970 



7?hii = 3.435 x 10" 30 -; ' U ^' H \_ 3 720 n c n H n (A24) 



Am \ - 376 



VHeii = k b T ai cU n e nH e u 



V i eIU = 2.748 x 10" 29 



rp 1.970 



(A25) 
(A26) 



[1-0+ (W . 
• Bremsstrahlung Cooling (Cen, 1992) [ergs cm" 3 s" 1 ] 
= 1.42 x 10" 27 5// \/Tl(nHii + riHeii + 4n H eiii)n e (A27) 
where gff = 1.5 is the Gaunt factor. 



• Compton Heating/Cooling (Haiman et al., 1996) 

[ergs cm" 3 s" 1 ] 

(A28) 



X = 1.017 x 10~ 37 T 7 4 (To - T 7 ) n c 



where T 7 = Tskgnd/K is the unitless background radia- 
tion temperature. 
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